A Hybrid Monte Carlo Method for Surface Growth Simulations 
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We introduce an algorithm for treating growth on surfaces which combines important features of 
continuum methods (such as the level-set method) and Kinetic Monte Carlo (KMC) simulations. 
We treat the motion of adatoms in continuum theory, but attach them to islands one atom at a 
time. The technique is borrowed from the Dielectric Breakdown Model. Our method allows us to 
give a realistic account of fluctuations in island shape, which is lacking in deterministic continuum 
treatments and which is an important physical effect. Our method should be most important for 
problems close to equilibrium where KMC becomes impractically slow. 

PACS numbers: PACS numbers: 68.55.Jk, 68.35.Fx, 81.15.Aa 
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Epitaxial growth on surfaces is of central importance 
both for applications and as a very interesting example of 
statistical processes out of equilibrium. We may idealize 
the process as the introduction of new atoms (adatoms) 
onto a crystal surface with flux, F; the adatoms then 
diffuse, with diffusion coefficient, D, nucleate islands or 
attach to existing islands. During the early stages of 
growth, the submonolayer case, the island size and shape 
distribution is a matter of substantial practical and the- 
oretical interest. 

The island growth process is commonly modeled by 
kinetic Monte Carlo (KMC) or continuum models. In 
KMC, internal noise processes are automatically repre- 
sented within the model and each adatom is represented 
individually. Therefore, when there are many adatoms 
(e.g. close to equilbrium) such simulations slow down 
considerably. A deterministic continuum model which 
represents the adatoms as a continuous fluid does not 
have this problem, and should be much faster. There has 
been considerable work in the development of such mod- 
els for epitaxial growth (see 0,0] and references therein). 
In some cases they have been quite successful, but some- 
times they do not reproduce experimental results. One 
reason for such problems is that deterministic continuum 
models neglect important fluctuations. In this paper we 
present a method of dealing with some fluctuations with- 
out giving up the advantages of a continuum treatment. 
We call this approach Hybrid Monte Carlo (HMC). The 
most important use of this method will be in cases where 
fluctuations are important, but which would be difficult 
to treat with KMC because of the presence of a large 
number of adatoms. A computation in a similar spirit 
has been given in • However, the present approach 



differs in a number of important ways. 

There are several sources of fluctuations in the growth 
process. The one we consider here is the fact that when 
atoms attach to islands they do so one at a time - that 
is, there is shot noise in the island growth process. This 
is important because island growth limited by diffusion is 
intrinsically unstable. The surface of the island will grow 
fingers due to the analogue of the well- known Mullins- 
Sekerka instability of metallurgy jj]. In the context of 
thin film growth the instability was discussed in detail by 
Bales and Zangwill The reason for unstable growth 
is easy to see: if a finger on the edge of an island starts 
to grow it will project out onto the terrace and be fed by 
more adatoms than the portions behind. The finger will 
grow longer, and be fed by still more flux, etc. Edge dif- 
fusion and other restructuring processes smooth out the 
fingers, and the final shape depends on the competition 
between unstable growth and smoothing. 

An extreme case of the diffusively unstable growth is 
represented by the Diffusion-Limited Aggregation (DLA) 
model of Witten and Sander Q. In this model all 
smoothing processes are neglected and growth takes place 
so slowly that one random walking adatom at a time is 
considered. DLA clusters are sprawling fractal objects 
with many branches which resemble some cases of island 
growth [alii- There is a variation of the model called the 
Dielectric Breakdown Model (DBM) 10] in which ran- 
dom walkers are not used. Instead the Laplace equation 
is solved outside the aggregate for a field, p, which repre- 
sents the probability density of walkers, and the growth 
algorithm is to add one particle at a time with probabil- 
ity proportional to dp/dn at the surface. Thus the DLA 
limit can be successfully treated by a model in which 
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the adatoms are a continuum. However, methods such 
as the level-set method cannot go to this limit and 
thus cannot produce dendritic islands which are seen in 
experiment. 

The unstable modes for the interface of the islands are 
present in level-set models, of course. However, the rea- 
son why the DBM limit is not achieved is that the growth 
process is represented by deterministically advancing a 
continuum interface according to the flux of the adatom 
fluid into the surface (see Eq. J3J below). In this algo- 
rithm the amplitude of the perturbations to a smooth 
interface are not correctly represented: they are given 
either by the initial conditions or by computer roundoff 
errors. In the experiment, however, there is a mechanism 
for feeding the instability: each adatom attaches not as 
a spread-out advance of the interface, but as an atom. In 
the case of DLA the result is that noise in the shape is 
present at all scales 01 an d does not average out. 

We should mention that there are other fluctuations 
which we will not treat. For example, there are density 
fluctuations in the adatom fluid which are important in 
the nucleation of new islands. There is a method 0] to 
treat this within level-set theory that we could employ. 
We will not consider such processes in this work, but 
rather look at the shape fluctuations of existing islands. 

The HMC algorithm goes as follows: we treat the is- 
lands on the surface as crystals containing discrete atoms 
which occupy the sites of a lattice. To illustrate the 
method we use a square lattice here. On the other hand, 
the adatoms are treated as a continuum whose surface 
density is p, and which is governed by: 



dtp 
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(1) 



In practice we solve this equation numerically on a dis- 
crete square grid which is commensurate with the crystal. 

Eq. (JJJ is solved with periodic boundary conditions 
on the edge of the system. The important physics of 
growth in incorporated into the boundary condition at 
the surface of the island. We put: 



dp 

D-^ = Tk±(p- Po) 



(2) 



Here k- and k+ arc the respective attachment rates on 
the upper and lower edge of the island boundary. 

In the case of irreversible growth only the first term in 
brackets in Eq. J5J) would be present. The other term, 
p , accounts for the detachment of atoms from the is- 
land, and depends on the position on the island bound- 
ary. Physically, the boundary condition must allow for 
faster detachment at corners, say, than at flat surfaces. 
We represent this in a way that allows us to compare 
directly with the bond-counting version of KMC: 



Po = exp(-n£ , /fc B T) 



(3) 



Here n is the number of nearest neighbor bonds that must 
be broken to completely detach (see below) the atom in 



question and add it to the adatom sea. The value of n 
depends on the environment of the detaching atom. We 
imagine that atoms can break bonds by moving along 
lattice directions. Also, E is the bond energy, and T 
the temperature. Of course, we can easily incorporate 
bonding to more distant neighbors. 

In a pure continuum model, the velocity of the island 
growth would be determined by mass conservation: 



v n — a D 



dp 
dn 



(4) 



where a is the lattice constant and [•] denotes the jump 
across the island boundary. Note that we can interpret 
Eq. J5J) in a way which is familiar in studies of crystal 
growth |F3j . By combining Eqs. (121 111 we find an equation 
for the density at the surface: 



P - P = {po - p) + av n 



(5) 



where a = l/(a 2 [fc + + &_]), and p is the equilibrium den- 
sity near a flat surface. That is, we are including both 
local equilibrium and kinetic terms in our boundary con- 
dition. The difference p — ~p is a measure of the number 
of dangling bonds on the surface, and thus of the curva- 
ture (upon coarse-graining). That is, the first term in Eq. 
(JSJ is related to the familiar Gibbs-Thompson boundary 
condition of crystal growth, and the second is a kinetic 
term. 

In HMC we implement Eq. Q in a way that includes 
fluctuations. Consider first a case where attachment is 
the only important process. Then we solve Eq. Q[2J) and 
compute the total flux onto the island boundary using 
Eq. |0 0} • When the total flux exceeds one atom then 
an adatom is attached to the boundary at random with 
the probability proportional to v n (exactly as in the DBM 
model.) In the case where detachment is also present we 
consider the surface to be partitioned into the part where 
the net flux is inward (growth), and attach atoms with 
probability density oc v n , and outward (detachment) and 
remove island atoms with probability density oc — v n . 

We have implemented HMC and compared the results 
to a KMC code on the same square lattice using near- 
est neighbor bonding. We use the hopping rate of an 
adatom to set the unit of time, and the lattice constant 
to be unity. We have two independent parameters namely 
D/F. and e = E/ksT We can also add edge diffusion, 
but here we have not done so. The KMC code is written 
using the method described in [J^j which takes advantage 
of the fact that there are only a few independent jump 
probabilities for a bond-counting model. We find that, as 
expected, at large e, the KMC code is much faster. How- 
ever, for e = 1.5 and with about 300 adatoms in a 40x40 
system the speeds are comparable. We note that in situ- 
ations such as heteroepitaxy |l5l | there are a great many 
independent probabilities to jump and KMC is slower. 

The interpretation of the number of bonds, n, is a bit 
delicate. In Fig. JJ| we show some examples of what we 



FIG. 1: a.) Detaching from a [10] surface can be done in 
one step, and breaks 3 bonds. b.)Fully detaching from a [11] 
surface in two steps also breaks 3 bonds. » 




mean by 'complete detachment'. For our square lattice, 
to detach an atom on a [10] surface we need to break 3 
bonds. However, for a [11] surface, in order to detach 
from the surface, an atom must first break two bonds, 
and then subsequently, one more - see Fig. Thus, 
for a [11] surface we set n — 3 as well because this cor- 
responds to the product of the probabilities of the two 
processes. In effect, we have coarse-grained. For other 
possible surface environments it is not difficult to tabu- 
late the correct value of n. 

This procedure may seem arbitrary, but we can justify 
it by its results. To this end, we show that HMC gives the 
correct shape for an island in equilibrium with adatoms. 
This shape is known exactly |l6( from a mapping to the 
2d Ising model. We did simulations with F = starting 
with a square island, and ran our code for a long enough 
time that the system seemed to be in equilibrium. In 
Fig. J2J we show some results superimposed on the exact 
results for various temperatures. The simulation results 
are ensemble averages; that is, we did 20 independent 
simulations and averaged the density of the island after 
shifting the center of mass of each one to be at the origin. 
The dotted line in Figure [3 is the contour line where the 
averaged island density is 1/2. 

We now present some more results to demonstrate the 
technique. If growth dominates detachment we should 
have DLA-like structures. Whether this occurs depends 
on T and D/F t l7.]. For low temperature, Fig. ©, we 
see the transition in a very clear way. Fig. shows a 
higher T case. Of course, all of these effects can be seen in 



FIG. 3: Irreversible (low T) growth of islands for two values 
of D/F. Top, D/F = 10 5 . Bottom, D/F = 10 2 . 
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FIG. 4: Island growth with D/F = 10 5 and e = 2.5 

KMC simulations. In other simulations (not shown), we 
have demonstrated that edge diffusion also smoothes out 
dendritic shapes, as expected. The virtue of our method 
will be to treat systems near equilibrium where the dy- 
namics of fluctuations are of interest and where there are 
a great number of adatoms, as in Fig. 

Another example is the thermal broadening of steps 
due to repeated attachment and detachment of adatoms. 
There is an extensive literature in this area. 0] , and the 
theoretical expectation [l8lll9| is that the thermal width, 
w, of a step should depend on the rate- limiting mecha- 
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FIG. 5: Thermal roughening of a pair of steps, average over 20 
simulations. Lower curve (dotted) is the average adatom den- 
sity, p. For early times, while the steps are moving, and p is 
changing, the scaling is w 2 oc t 2/l3 , kinetic roughening. Later 
there is thermal roughening with w 2 oc t 1 ^ 2 (evaporation- 
condensation kinetics) crossing over to terrace-diffusion ki- 
netics, w 2 oc t 1 / 3 . The straight lines, from left to right, have 
slopes 2/3, 1/2, and 1/3. 



nism for step motion. In our model, without surface dif- 
fusion, this will be the either detachment from the step or 
diffusion on the terrace. In these cases w 2 oc t 1 / 2 , i 1 / 3 , re- 
spectively. We show, in Fig. |J5}, w 2 (t) for a pair of steps 
one atomic layer high which cross a 100 x 100 terrace. 
We see indications of both of the expected behaviors at 
late times. The early time behavior is kinetic roughen- 
ing j^Jj, w 2 oc t 2 / 3 . This is because we started with no 
adatoms, and, initially, the steps were retreating. 

The HMC technique should be most useful in situa- 
tions where there is a large separation of time scales be- 
tween the diffusion of the adatoms and fluctuations of 
island boundaries. Other examples for which it might be 
used are in studies of homogeneous nucleation of large 
islands near equilibrium on surfaces \2l\ . Another exam- 
ple is in the catalysis of the reaction CO+O — > CO2 on 
a noble metal. The diffusion of CO is very fast in this 
case, and makes a realistic KMC simulation very difficult. 
Evans and collaborators have used an approximation 
where pco is constant where ever there is no O on the 
surface. Our approach (with suitable changes in bound- 
ary conditions to allow for the reaction) should give a 
better account of the kinetics. 

Finally, we should mention the case of heteroepitaxy. 
There have been a number of KMC simulations of this 
very important process [l5ll23ll2^ |. This is a very difficult 
problem to attack numerically. The time-consuming step 
in such calculations is not the dynamics of the adatoms, 
but rather the elasticity computation. Nevertheless, our 
approach to the adatoms offers a number of advantages, 
and we intend to pursue this in a future publication. 
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